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Abstract. As we approach solar convection simulations that seek to model the interaction 
of small-scale granulation and supergranulation and even larger scales of convection within the 
near-surface shear layer (NSSL), the treatment of the boundary conditions and minimization 
O ' of sub-grid scale diffusive processes become increasingly crucial. We here assess changes in 

the dynamics and the energy flux balance of the flows established in rotating spherical shell 
C/^ . segments that capture much of the NSSL with the Curved Spherical Segment (CSS) code using 



two different diffusion schemes. The CSS code is a new massively parallel modeling tool capable 
of simulating 3-D compressible MHD convection with a realistic solar stratification in rotating 
spherical shell segments. 



OO . 

, 1. Introduction 

' The solar differential rotation profile exhibits prominent radial shear layers near the top and 

PsJ ■ bottom of the convection zone [1]. The near-surface shear layer (NSSL) occupies the upper 5% 

of the Sun by radius, whereas the tachocline begins near the base of the convection zone. The 
dynamics of the NSSL are governed largely by vigorous granular-scale convection that is driven 
by radiative cooling and large thermodynamic gradients. The collective interaction of these 
granular-scale flows (average sizes of 1 Mm, lifetimes of 0.2 hr) is a major component in the 
^ ' formation of supergranular (15-35 Mm, 24 hr) and mesogranular (5-10 Mm, 5 hr) scales 

Given this wide range of spatial and temporal scales, we currently cannot simultaneously 
model hundreds to thousands of super granules, solar granulation, and deeper flows. However, 
we may still be able to characterize the influence of granulation and deep global flows on the 
NSSL by coupling CSS with global convection models from below 0] and with surface convection 
simulations from above [e.g.H;[^. Such a coupling, whether it is statistical or direct, requires 
a careful treatment of diffusion and boundary conditions [H]. For instance, low diffusion is 
necessary to preserve the spatial structure and advective timescales of the downdrafts flowing 
into the CSS domain from the surface convection above. Given that our grid is five times coarser 
than that of typical surface convection simulations, this is no easy task. 

We have conducted two numerical simulations in a 20° square patch centered on the equator 
that encompass most of the NSSL and some of the deep interior, rotating at the solar rate. 
These simulations explore the effects of turbulent-eddy and slope-limited diffusion schemes with 
an open lower radial boundary and closed upper boundary; as such they are identical except for 
the diffusion scheme. The governing equations and numerical approach used in solving them are 
briefly discussed in ^ The turbulent-eddy and slope-limited diffusion schemes are detailed in 
J}3l The dynamics of the fiows established in these simulations are examined in ^21 




Figure 1. Radial velocities {ur) at 0.99 Rq in (a) the slope-limited diffusion Case 2, and (b) 
the turbulent-eddy diffusion Case 1, and at 0.95 Rq in (c) Case 2 and (d) Case 1. Dark tones 
denote downflows and light tones denote upflows. Scaling values are -1.34 kms-\ 0.59 kms"^ 
for (a) and (b) and -0.73 kms~^, 0.24 kms~^ for (c) and (d), where positive Ur points outward. 



2. Formulating the Problem 

The spherical segment domains used in our simulations involve large portions of the Sun's 
inherent spherical geometry, which is necessary to properly capture the effects of rotation on 
supergranular scales. We use the Curved Shell Segment (CSS) code to model the 3-D dynamics 
below the solar photosphere. CSS is a modeling tool that solves the compressible Navier-Stokes 
equations in rotating spherical segments (sj. To simulate the larger scales of motion that are 
likely to occur in the solar convection zone, a large-eddy simulation (LES) model is employed. 
The scales that are not explicitly computed in these simulations are parametrized and included 
in a sub-grid scale model of turbulent transport and diffusion. The governing equations solved 
in CSS are: 
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The symbols p, u, P, T, S, V, e, 5 are the density, velocity, pressure, temperature, specific 
entropy, viscous stress tensor, stress tensor, and identity tensor respectively; <I> is the viscous 
heating term; the entropy perturbations are S' = S — {S); the equation of state is P = p'^e^^'^^; 
fl is the angular velocity of the rotating frame; p and ks are the turbulent-eddy diffusion of 
momentum and entropy, and Kr is the radiative thermal diffusion; g is the local acceleration due 
to gravity, where p, ks, Kr, kq, and g are functions of radius only. 

The governing equations (Equations 1-3) are evolved on a uniform spatial mesh. Temporal 
discretization is accomplished using an explicit fourth-order accurate Runge-Kutta time-stepping 
scheme. The spatial derivatives are computed using a modified sixth-order compact finite 
difference scheme. A 3-D domain decomposition divides the full spatial mesh into sub-domains. 
The boundary information necessary to compute spatial derivatives is passed between nearest- 
neighbor sub-domains using MPI, while computations are shared among the master and slave 
cores within a supercomputer node using OpenMP. 

Currently there are two diffusion schemes implemented in CSS, a turbulent-eddy diffusion 
(TED) (Equation 3-4) and the slope-limited diffusion (SLD) described in (g3]). In the TED 
scheme, the momentum and entropy diffusivities {p and ks) are calculated based upon the 



desired Rayleigh number at the upper boundary with the constraint that the Prandtl number 
and the dynamic viscosity be constant throughout the domain. There are no directly comparable 
diffusive parameters in the SLD scheme as the eddy diffusion terms have been dropped. One 
may still control the level of diffusion, however, with specific choices of which slope-limiter is 
employed and what characteristic velocity is used at cell interfaces. Figure 1 depicts snapshots 
of typical flows at two radii from the two diffusion schemes. 

3. Slope-Limited Diffusion In CSS 

Slope-limited diffusion has many possible formulations, and we have chosen to use a scheme 
similar to that found in |^. SLD is based upon a piecewise linear reconstruction within a 
finite volume (cell) centered on each gridpoint of the solution at each time step. The linear 
reconstruction leads to a solution that has discontinuities at the cell edges. SLD essentially 
acts to minimize these discontinuities so that the numerical scheme remains stable. While this 
model of diffusion is not as physically motivated as the eddy diffusion model, it holds many 
computational advantages. Indeed, with SLD in CSS, a higher level of turbulence is achieved 
for a given resolution and there is no need to evaluate second-order derivatives, which decreases 
memory usage by a factor of three and the execution time per time step by a factor of two. 

The slope of the linear reconstruction of the solution is given by a ratio of the downwind 
and upwind cell-center differences (r^ = AiUj/Ariij). This slope is "limited" by a function 0(rj) 
that belongs to a class of slope limiters that yield total variation diminishing solutions As in 
we use a linear combination of two such slope limiters, the minmod and superbee limiters. 
Using the reconstructed slope, values of the primitive variables are computed at the cell edges 
as uf = Ui ± ^(f> (vi) ArUi. The diffusive flux at a cell edge is /f = ^cfgfl3f6uf, where gf is 
a geometric factor and c| is a characteristic velocity at a cell edge, and (5n| is the difference 
of the left and right reconstructed values at a cell-edge. To avoid artificial steepening, a local 
diffusion coefficient f5f is constructed such that when SufA^Ui > 0, f3'^ = {5ul/ AeUi)^ and is zero 
otherwise. This preserves steep gradients and yields a fourth-order upwind diffusion scheme [H]- 

With the diffusive fluxes at each edge of a given cell as above, the diffusion at the center of 
a cell is Di = i^fl — ff) / Axi, where Axi is the grid spacing. This is added to the solution after 
the full Runge-Kutta time step as Ui = Ui + AtDi. Within our simulations, using the true sound 
speed in the proves to be overly diffusive. So, for the computation of the cf , the sound speed 
throughout the computational domain is flxed to a fraction (~ 1/10) of the surface sound speed. 

4. Comparing Two Diffusion Schemes 

The simulations encompassing a 20° x 20° domain are relatively low resolution with 128 radial 
mesh points and 256 points in latitude and longitude. A stellar evolution code is used to establish 
a realistic initial stratification for the simulations. Since a perfect gas is assumed, the He and 
H ionization zones cannot currently be simulated. Thus, the impenetrable upper boundary is 
taken to be 0.995 Rq in order to exclude most of the radial extent of these zones. The permeable 
lower boundary is placed at 0.915 Rq, yielding a density contrast across the domain of 400. Each 
case has the lowest diffusivities that allow numerical stability in the TED [Case 1) and SLD 
{Case 2) regimes. Using a constant with depth in the TED case is rather restrictive in that 
the diffusion is larger than necessary away from the upper boundary, but for this preliminary 
analysis it is sufficient to characterize the differences between the two schemes. 

Turning to the radial velocity patterns sampled in Figure 1, there is a striking shift toward 
smaller spatial scales in Case 2 with considerably less power at large scales. At 0.99 Rq the 
peak of the spatial power spectrum in Case 1 is around 40 Mm. In Case 2, however, there is a 
broad peak at scales around 20 Mm, more closely matching the supergranular scales of the Sun. 
As one descends deeper into the simulations (Figures lc,d), the flows appear significantly more 
turbulent in Case 2 than in Case 1. 
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Figure 2. (a) PDFs of radial velocities at 0.99 Rq (solid line) and 0.95 Rq (dashed), for Case 
2 (black) and Case 1 (red). Time-averaged radial energy flux balances for Case 1 (b) and Case 
2 (c), with total flux in purple, radiative flux (blue), enthalpy flux (green), kinetic energy flux 
(orange), and acoustic flux (red). 

Indeed, the downflowing structures of Case 2 become much more narrow at depth with many 
more isolated plumes. The greater turbulence also manifests itself in larger extrema of each of 
the physical variables. For example, the radial velocities at 0.95 Rq in Case 2 have extrema 
that are twice that of Case 1 (Figure 2a). In Case 1, a precipitous drop in the upflowing radial 
velocity probability density distribution (PDF) outside of about 0.1 kms~^ is indicative of just 
how uniform the upflows are. In stark contrast, the radial velocity PDF of Case 2 has twice 
the dynamic range, emphasizing that the upflows are indeed much more structured than those 
of Case 1. While the radial mass flux PDF has higher velocity wings that contribute only one 
percent to the overall distribution, the bimodal nature of the temperature perturbation PDF 
has the effect of significantly amplifying the importance of these outlying values. Thus, the net 
effect of the larger wings in these PDFs is to shift the mean of the enthalpy flux in Case 2 to 
a higher value. The additional velocity in both the horizontal and radial flows greatly alter the 
kinetic energy distribution, as in Case 2 where the mean kinetic energy has nearly doubled and 
correspondingly increased the magnitude of the kinetic energy flux (Figure 2c). 

A tentative next step is to have both diffusion schemes active at the same time. This will 
allow a systematic study on the convergence of solutions as the TED coefficients are lowered to 
solar values and the SLD takes over as the primary diffusion. 
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